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We study the effects of motor-generated stresses in disordered three dimensional fiber networks 
using a combination of a mean-field, effective medium theory, scaling analysis and a computational 
model. We find that motor activity controls the elasticity in an anomalous fashion close to the 
point of marginal stability by coupling to critical network fluctuations. We also show that motor 
stresses can stabilize initially floppy networks, extending the range of critical behavior to a broad 
regime of network connectivities below the marginal point. Away from this regime, or at high stress, 
motors give rise to a linear increase in stiffness with stress. Finally, we demonstrate that our results 
are captured by a simple, constitutive scaling relation highlighting the important role of non-affine 
strain fluctuations as a susceptibility to motor stress. 
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The mechanical properties of cells are regulated in 
part by internal stresses generated actively by molecu- 
lar motors in the cytoskeletal filamentous actin network 
[l|. On a larger scale, collective motor activity allows 
the cell to contract the surrounding extracellular ma- 
trix, consisting also of biopolymer networks. Experi- 
ments show that such active contractility dramatically 
affects network elasticity, both in reconstituted intracel- 
lular F-actin networks with myosin motors J2h5| and in 
extracellular matrices with contractile cells [f|. The dy- 
namics and elasticity of active biopolymer networks have 
been studied theoretically using long-wavelength hydro- 
dynamic approaches as well as affine models 043 ■ These 
approaches, however, fail to describe highly disordered 
networks. There is also experimental evidence that cy- 
toskeletal networks may be unstable or only marginally 
stable in the absence of motor activity In such 

cases, networks are expected to be governed by highly 
nonuniform, soft or floppy modes of deformation that 
may lead to a fundamental breakdown or failure of con- 
tinuum elasticity (llj . Importantly, motor- induced con- 
tractile stresses can be expected to couple to these soft 
modes, giving rise to a nonlinear elastic response that is 
distinct from the nonlinearities arising from single fiber 
elasticity that have been considered in previous models. 
Moreover, such a coupling to local soft modes of the net- 
work may call into question the equivalence of internal 
(motor) and external stress, a tacit assumption in the 
analysis of recent in vitro experiments [H, 0| • 

Here, we introduce a simple model to study the effects 
of motor generated stresses in disordered fiber networks. 
Networks are formed by crosslinked straight fibers with 
linear stretching and bending elasticity. These fibers are 
organized on a face centered cubic (FCC) lattice in which 
a certain fraction of the the bonds can randomly be re- 
moved. Motor activity is introduced by contractile force 
dipoles acting between neighboring network nodes. We 
find that motors can stabilize the elastic response of oth- 
erwise floppy, unstable networks. The motor stress also 



controls the mechanics of stable networks above a char- 
acteristic threshold, in the vicinity of which the network 
exhibits critical strain fluctuations. We develop a quan- 
titative effective medium theory to describe the elastic 
response of these systems. Interestingly, the network's 
stiffness is controlled by a coupling of the motor induced 
stresses to the strain fluctuations. This coupling gives 
rise to anomalous regimes at the stability thresholds, at 
which network criticality is reflected in both divergent 
strain fluctuations and anomalous dependences of the 
network mechanics on stress. In these critical regimes, 
the shear modulus depends nonlinearlyon both motor 
stress and single filament elasticity (fl fl^4T3 |. Interest- 
ingly, this dependence on internal motor stress differs 
qualitatively from that of an applied external stress. 

A key parameter that characterizes fiber networks is 
the mean coordination number, z. Although the net- 
work is connected above a threshold z = z C ond — 2, 
it only becomes rigid above a higher rigidity threshold 
Zb ~ 3.4. This threshold is due to the bending rigidity 
of the individual fibers and it lies below the central-force 
(CF) rigidity threshold, z c f ~ 6, for a spring-only net- 
work. In general, when some fraction of the bonds are 
under stress, additional constraints are introduced [l5) . 
More formally, these constraints appear as scalar terms in 
the Hamiltonian [l(| . These additional stress-constraints 
may shift the various rigidity thresholds in the system. 
In random spring networks, for example, this can be re- 
alized by applying finite network deformations; this has 
been studied in spring networks 

EMU 

where the ac- 
tual rigidity threshold shifts continuously to lower values 
with the applied external strain. Under such external de- 
formations, the internal stress is free to adopt the most 
favorable distribution. By contrast however, motors im- 
pose a fixed distribution of internal stress, which may 
lead to a qualitatively different network mechanics. 

To provide insight into the elasticity of fibrous net- 
works with contractile internal stresses, we use a model 
of fibers organized on a FCC lattice. By removing lattice- 
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bonds with a probability 1 — p, we tune the average co- 
ordination number, z = Zp, where Z = 12 for the undi- 
luted lattice. Motors are introduced as contractile force 
dipoles and are inserted randomly with a probability q. 
The fibers are modeled as linear elastic beams with a 
stretching modulus /i and bending rigidity K. Using units 
in which £q = /i = 1, the total energy can be written as 
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r,' and denotes the position of z'th 



node and Ptj = 1 for present bonds or 
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moved bonds. The first sum extends over neighboring 
pairs of vertices. The crosslinks themselves do not con- 
tribute a torsional stiffness and, thus, the second sum 
only extends over coaxial nearest neighbor bonds on the 
same fiber. The last term represents the work performed 
by the motors, where Qij = 1 if a motor acts between 
nodes i and j and Qij = otherwise. 

To develop a mean-field, effective medium theory 
(EMT) that captures the disordered nature of this 
model — including internal stresses — we extend the the- 
ory for the linear mechanical response of disordered 
spring networks (20l - f2^ | . In our EMT approach we ignore 
the bending contribution (k = 0), allowing us to circum- 
vent the difficulties involved in an EMT with three-point 
bending interactions [ill, HI HI- Our EMT is based 
on a mapping between the disordered and an ordered 
network with an effective elastic constant, yet with the 
same underlying lattice geometry and under the same 
internal stress as the original disordered system. The ef- 
fective elastic constant, /I (cm), is determined by a self- 
consistency condition; the local distortion in the effective 
medium induced by replacing a bond, selected randomly 
from the disordered system, should vanish on average. 
For a general disordered network this procedure yields 
an implicit expression for the effective stretch modulus 
(see Appendix) 



fJ'EM ~\~ fJ>ij (Ji (cm) 



V( f x ij )dfi ij =0, (2) 



where Hem is the displacement of a bond in the unper- 
turbed effective medium due to a unit force acting along 
the bond, \Xij is the stretching modulus between nodes i 
and j and V (fiij ) is the probability density of the moduli 
in the disordered system. For the case of a diluted lattice 
considered here, V (fJ>ij) = pS (fXij — 1) + (1 — p) S (/%), 
we find the EMT shear modulus 
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FIG. 1: (a) Shear modulus is shown vs. motor stress for 
different values of the coordination number z. Markers repre- 
sent the numerical data (see legend in (c)). Thin, solid lines 
represent the results of the EMT for the same values of z 
while upper lines correspond to higher values of z. Lines with 
numbers represent the corresponding slopes and thick grey 
area contains the afRne predictions for all presented coordi- 
nation numbers, (b) G is shown vs. z for different values of 
<jj\/. Markers represent numerical data (see legend). Thin, 
solid lines represent the results of the EMT for the same val- 
ues of om while upper lines correspond to higher values of 
om- Thick grey area contains the affine predictions for all 
presented values of <tm- (c) The differential non-affinity pa- 
rameter is presented vs. motor stress for different values of z 
(see legend), (d) 8F is presented vs. z for different values of 
am (see legend in (b)). For all data sets k = 10~°. 



where, within the EMT, om = V8qf. 

While the full expression for /I((Tm) is long (see Ap- 
pendix), the scaling predictions of the EMT are simple. 
Even below the central- force isostatic point, z c f, motor 
activity induces a finite shear modulus. Far from z c f, 
G ~ Go + cjvf, where Go is the shear modulus of the 
unstressed network 1 . By contrast, close to z c f there is an 
anomalous scaling regime G ~ er^ 2 /! 1 / 2 . 

To test the implications of the EMT, we perform sim- 
ulations of fiber networks with finite bending rigidities. 
The shear modulus, G, is determined numerically by 
applying a shear strain along the Ill-plane using Lees- 
Edwards periodic boundary conditions and energy min- 
imizations are performed by a conjugate gradient algo- 
rithm . First, we consider the high motor density 
limit q ~ 1. The EMT prediction is in good quantitative 
agreement with the numerical results over a broad range 
of network connectivity and motor stress, as shown in 
Fig. [1] a,b. Since we neglected the contribution of fiber 



1 In all scaling relationships with additive contributions, unknown 
numerical prefactors are omitted. 
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bending energies in the EMT, it fails in the regime where 
G is governed by k. In addition, in the vicinity of z c f for 

cm k, we find a mixed regime, G ~ ff^~ v <r v M , where 
y' = 0.4, whereas in the EMT y' = 0.5 (Fig. EJ). This 
mixed regime is similar in nature to the n-(i coupled me- 
chanical regime around z c f in unstressed fibrous networks 
[Tlj . More generally, such coupled regimes arise in the 
vicinity of a stability threshold, when there are additional 
interactions or fields that stabilize the network below the 
threshold Thus, in this model the motor stress acts 
as an external field. In fact, as may be expected, another 
anomalous regime is observed in the simulations at the 
bending rigidity threshold, G ~ fv~~ y cr v M , with y = 0.6 
(Figs. [Hand©. 

We gain additional physical insight into the elastic 
properties of active networks with a scaling argument 
we estimate the amount of work that is performed by the 
motors when the system is sheared. The characteristic 
deformation of a single bond in such a network will be 
such that it avoids energetically costly stretching contri- 
butions. Such deformations are oriented perpendicularly 
to the direction of the bond: the nonaffine contribution 
to this deformation can be estimated by 5u± ~ "fs/ST, 
where the differential nonaffinity parameter is defined as, 



sr = ±- 2 ((su k -6ufy 
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Here Suk is the displacement of node k under an infinites- 
imal external shear 7, Su^f is the affine prediction and 
the average is taken over all network nodes. Interest- 
ingly however, this is not the only relevant contribution 
to the deformation of the bond. The component of the 
affine deformation perpendicular to the bond does not 
contribute to bond-stretching energies to harmonic order 
and, thus, is not avoided. Importantly however, this de- 
formation does contribute to the motor work. Therefore, 
the total work performed by the internal stress resulting 
from such deformations scales as 8W ~ aMl 2 ST + aMj 2 , 
implying the following relationship for the shear modu- 
lus, 



G ~ Go + (jmST + <jm. 



The non affinity parameter, ST(ctm, z, k), depends on 
the system's parameters as shown in Fig. [T}j,d. To con- 
firm the prediction of Eq. ([5]) we plot G — Go — §0m vs. 
Cm ST and find that all data collapses on to the same 
curve with a linear dependence, as shown in Fig. [3Ja). 
Interestingly, the scaling prediction in Eq. (|5j) suggests 
that ST can be interpreted as a susceptibility of the shear 
modulus to the internally generated stress. Moreover, ST 
shows a strong increase close to both rigidity thresholds 
(Fig. Q] d) , implying a large susceptibility to om when 
the system is marginally stable. However, at these sta- 
bility thresholds ST acquires a strong dependence on <tm . 
This can be understood by considering ou as an external 
field that restores rigidity and suppresses the divergence 
of the strain fluctuations ST ~ |Az|~ (for <jm = 0); this 
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FIG. 2: Demonstration of the anomalous regimes, (a) Shear 
modulus separated to the linear part, 5/6<tm (dashed line), 
and the rest, G — 5/6<tm (circles) is shown for 2 = 3.35 (close 
to the bending rigidity percolation point, Zb — 3.36) and k = 
10 -2 . In the small stress limit the G — 5/6<tm part is found 
to scale as a]f y with y ~ 0.6 (solid line), (b) The same 
analysis as in (a) is shown for 2 = 5.9 (close to the central- 
force rigidity percolation point, z c f — 5.64) and ft = 0. In 
the small stress limit the G — 5/ 6a m part is found to scale as 
o~\j y with y' ~ 0.4 (solid line), in contrast to the mean-field 
prediction (indicated by the dashed-dotted line) that scales 
as y/aM in the small stress limit, (c) The value of ST (circles) 
is presented for the same set of parameters as in (a). In the 
small stress limit it scales as a~^f (solid line), (d) The value 
of ST (circles) is presented for the same set of parameters as 
in (b). In the small stress limit it scales as cr^/ (solid line). 



implies a dependence of the form ST(zb) ~ {k/ctm) v and 
ST(zcf) ~ ((i/ctm)~ v and, taken together with Eq. (0, 
explains the origins of the anomalous regimes, where 



G(z b ) 



G(z ci ) ~ vfi 
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We verified the internal consistency by determining the 
exponents y and y' at the two stability thresholds from 
both the scaling of G and ST with <tm, as shown in Fig. 



The schematic phase diagram for the high-motor den- 
sity limit is shown in Fig. [3(b). Away from the stabil- 
ity thresholds the shear modulus scales linearly with the 
active stress. This is in contrast with the stiffening be- 
havior of externally deformed networks, for which the de- 
pendence of the differential elastic modulus goes as the 
square root of the external stress [l8l . [l9j . Thus, there 
is not necessarily a quantitative correspondence between 
internally and externally stressed networks, in contrast 
to suggestions in prior work 

Finally, we explore the role of inhomogeneity in the 
distribution of active motors, which shows that critical 
behavior is not limited to the critical points associated 
with rigidity percolation. We model inhomogeneous mo- 
tors by considering the range q < 1 for different values 
of z well below the rigidity percolation point, z < Zb- 
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FIG. 3: (a) Collapse on the data presented on Fig. QJa) based 
on Eq. (0. Red line represents linear dependence, (b) The 
schematic phase diagram for the rigidity of random spring 
networks under an internal stress ctm- 



1(T 
G 

10" 6 

1Cf 7 
10° 

sr 



□ 2.0 
♦ 2.4 

O 3.0 



(c) 



♦ g 



10 1<T 

cm 



10" 

cm 
10- 5 

10" 6 
10" 7 

0.6 

sr 

0.4 



(b) 



o° > D 
I 5 , 



Oo 



(d) 



In this case the motors only induce a macroscopic stress 
when the motor density exceeds a z-dependent threshold, 
q c (z), as shown in Fig. Hfb). Concurrent with the de- 
velopment of a macroscopic stress, the network acquires 
a finite shear rigidity. Near the threshold q c . the motor- 
induced stress falls significantly below the mean-field pre- 
diction ((7m = V&qf) and depends non-linearly on q. In- 
terestingly, in this regime (a m *C V8q c f) the nonaffine 
fluctuations become large (see Fig. Bid)), diverging with 
motor stress with an exponent close to —0.2, as shown 
in Fig. [He). Such a divergence, taken together with Eq. 
([5|) implies an anomalous, sub-linear scaling of the shear 
modulus with the motors stress with the exponent 0.8. 
Indeed, as shown in Fig. HJa), the stiffening of the shear 
modulus clearly deviates from the mean-field predictions 
and scales sublinearly with the motor stress with an 0.8 
exponent, consistent with Eq. ([S]), even when the mean 
coordination number of the network is well below the 
rigidity percolation point. Thus, criticality in the form 
of a divergent susceptibility is characteristic of floppy sys- 
tems below the rigidity percolation point. 

This work demonstrates that motor activity controls 
the elastic properties of disordered networks by coupling 
to the differential non-affine fluctuations in the deforma- 
tion field. This coupling makes elastic deformations more 
affine and stabilizes the network. Far from the elastic 
critical points this coupling leads to linear stiffening as a 
function of the motors stress, as has been observed in sev- 
eral studies of prestressed elastic networks [27| ■ How- 
ever, close to the elastic critical points, where the non- 
affine fluctuations diverge, this coupling leads to anoma- 
lous regimes, where the shear modulus scales sub- linearly 
with the motors stress. Similar stress-stiffening of floppy 
networks below marginal stability is also found beyond 
a threshold in the motor density, indicating that a sur- 
prising generality of critical fluctuations and divergent 
susceptibility for systems below the usual rigidity perco- 
lation point. 



FIG. 4: The role of inhomogeneity in the motor distribution 
for different values of z (see legend on panel (a) for every plot 
in this figure), (a) Shear modulus is presented vs. the nor- 
mal stress induced by the motors. Solid line represents the 
power law of 0.8 while the dashed line contain the mean-filed 
predictions, (b) The normal stress vs. the motor probability 
occupation. Solid line represents the mean-field prediction, 
om = V8qf. (c) The differential non-affine fluctuations mea- 
sure for the same data as in (a) is presented vs. ctm. The 
solid line represents the power law of —0.2. (d) SF vs. q. For 
all the data in this figure / = I0~ 4 and k = I0~ 5 . 
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Appendix A: Mean-field approach 

The nonlinear EM approach developed here is based 
on a scheme to construct a mapping from the internally 
stressed lattice network with disordered spring constant, 
/ijj, with probability density P(fiij), onto a perfect lat- 
tice system with uniform bond stiffness with the same 
stress, om- The stress of the network appears due to 
motors' force / that acts between each pair of neighbor- 
ing nodes of the network. This mapping is realized by 
an effective uniform central force interaction, /i.y — > £t. 
The effective parameter, J1(ctm), is determined by a self- 
consistency requirement: replacing a random bond in 
the uniform EM under stress with a bond drawn from 
the original probability density, P (/ijj ) , results in a local 
fluctuation in the deformation field, which vanishes when 
averaged. In addition, we assume that the fluctuations 
of the deformations are small compared to the distance 
between crosslinks. This approach leads to an integral 
equation, representing a disorder average (Eq. (|A15[0 . 
from which the effective parameter /I(ctm) can be deter- 
mined. In the following we ignore the contribution of the 
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bending stiffness of the network filaments. Due to this 
assumption the presented mean-field approach fails in the 
regime where the elasticity of the network is dominated 
by the bending stiffness. 

1. Effective medium theory 

We apply the EM theory method to a network sub- 
jected to a uniform internal compression resulting in a 
macroscopic isotropic stress <jm- Similarly to Ref. [20l |. 
we calculate the effective spring constant using the self- 
consistency requirement. 

The position of a crosslink (network's node) i is given 
by Ri = R? + Ui, where R° is the position in the un- 
stressed configuration and Uj is the displacement field. 
In order to calculate the elastic constant we apply an in- 
finitesimal external expansional/comressional strain e to 
the network. The affine displacement due to the applied 
strain is given by uf s — Uj S = ery , where ry is the vector 
from R? to in the undeformed reference state. Here 



The first term represents the expansion/compression en- 
ergy of the affine response, while the other terms cor- 
respond to the energy difference due to the non-affine 
deformation of the stretched/compressed bonds. 

The expansion of the whole network corresponds to the 
global constraint 

Yj*ij = 0. (A5) 
To investigate the elastic behavior of the model in Eq. 



where Ay = — (/ + e/Z) ry . To calculate the effective 
parameter jl we demand self-consistency of the EM [2(| ■ 
The self-consistency requirement in this context can be 
formulated as follows: the non-affine displacement in- 
duced by the replacement of a single bond in the EM 
vanishes on average, 



we allow for non-affine displacements 

v, = u, - uf , (Al) 
The Hamiltonian of the network is given by 

(ij) (ij) 

where Ry = Rj — R^ and the sums extend over neigh- 
boring pairs of vertices. We assume that the resulting 
non-affine relative displacements of neighbouring nodes 
i and j are much smaller than the distance between the 
nodes, 

|v„-|EE|v i -v J |«|R°-R J °| = l. (A3) 

Thus, we can expand the Hamiltonian around the affine 
strain configuration (small i>y). U p to second order in 
va and first order in e we arrive at [la, Il7fl 



(A4) 

I 

(|A4|) . we set up an effective medium theory. In the EM 
approach, we mimic the disordered system by the regu- 
lar one with an effective parameter, i.e. /iy — > Jl and the 
same stress cm- In other words, the EM network may 
be globally expanded by applying the force that assures 
mechanical equilibrium for the affine, vy = 0, configura- 
tion. Thus, the EM system has the Hamiltonian, given 
by 



(A6) 



Here, the average is taken over the distribution of the 
nm bond in the original disordered system, i.e. accord- 
ing to the probability density P (fx nm ). To calculate the 
displacement v„ m after the replacement we solve the per- 
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Hem = 5^/ + (/ + e(i) Vy • ry + -(£-/) (vy • ry) 2 + -/vy + Ay ■ vy 

I 



(v nm ) = 0. 



(A7) 
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turbed EM Hamiltonian that is given by 
Hem + ^ (f^nm — /•*) ( v 

(A8) 

In the configuration that minimizes the energy, the dis- 
placement of the nm bond is given by 



HEM + A*nm — H ' 



(A9) 



where hem is the displacement of the nm bond in the 
unperturbed EM network due to a unit force acting along 
the nm bond. 



a. The calculation of [i em 

In this Section we calculate hem — the displacement of 
the nm bond in the unperturbed EM network due to a 
unit force r nm acting on the nm bond. 

The dynamical matrix of the unperturbed EM Hamil- 
tonian (|A6[) is given by 



Di 



~(H~ f) Tij <g> Tij + 

£ Km - /) r u ® *ij 



fl i + j 
f/I] <=j 



(A10) 



where I is the unit tensor and (g> is the external product. 
The Fourier transform of D is given by 

D(k) = Y. l) :- , ' kr = 

ij 

= (0* - /) r (8 r + /I) (1 - e ikr ) (All) 

r 

where r runs over all unit bond vectors. The unit force 
acting on the nm bond is given by 

so that its Fourier transform is 

f (k) = J2^ lk B " = r nm (1 - e* k r -) . (A13) 



Thus the Fourier transform of the displacement field is 
given by 



v(k) = -i?- 1 (k).f(k) 



(A14) 



The displacement of the nm bond due to the unit force 
is 



f-EM 
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]Tv (k) (e" ik - r - - 1) = -J2*nm ■ f (k) D- 1 (k) (e- 
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For a highly coordinated lattice the sum over r may be sum over k is dominated by the small k • r <C 1 values, 
well approximated by the integral over the sphere that hem may be approximated by 
includes all the neighbouring crosslinks and, since the 
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dN (h - /) [ $ (k • r) 2 d^r (r <8 r + J^] 

_/ / l d-i 

d(u-f) [ -J- + J- + -J- + J— 



6. An effective elastic constant 



Given Eqs. (|A7IA9IA15|) . the self-consistency Eq. 
(|A7j) leads to the following equation for the effective pa- 
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rameter 2 



2. Diluted regular networks 



Hij - n(a M ) 

fJ-EM + Hij - M ((Tm) 



/';//,,)<///,, 0. (A15) 



The approach presented in this section allows one to 
calculate the elastic parameters of a system with a given 
topology and elastic constant distribution in the nonlin- 
ear elastic regime. Eq. (|A15|) may be solved numeri- 
cally for any realization of the spring constant probability 
density, P(iXij). Knowing the effective spring constant, 
ju, one obtains all the elastic constants of the network 
and the relation between the motors' applying force, / 
and the global normal network's stress, o~m- In the next 
section we demonstrate the presented method using the 
particular example of diluted regular networks when Eq. 
(IA151) can be solved analytically. 



In this Section we use the mean-field solution presented 
above using for particular example of bond-diluted reg- 
ular networks. The probability density for the spring 
constants for such a network is given by 

P ilHj) = Jf<5 - M ) + (l - |) 5 (/xy) . (A16) 

Networks of this kind are referred to as diluted spring 
networks or the central-force elastic percolation model. 
The linear elastic response of diluted, unstressed lattices 
has been extensively studied (20l . [28| . Here we show how 
these results generalize for internally stressed networks. 
In this case the Eq. (|A15|) becomes 
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z/Z 
1 - z/Z' 



(A17) 



The solution of this equation provides the spring constant 
of the effective medium, /I, and, therefore, knowing the 



geometry of the original lattice one can easily calculate 
all the elastic constants of the network. 
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